package statistics;

import java.math.BigDecimal;
import edu.rit.numeric.Statistics;

public class R_Contrasts extends AbstractStatisticTests implements IStatisticalTest 
{	
	
	static BigDecimal[] lnFactorialList = new BigDecimal[5000];
	static {
		lnFactorialList[0] = new BigDecimal("0");
		for(int i=1; i<5000; i++) {
			BigDecimal number = new BigDecimal( Math.log( i ) );
			lnFactorialList[i] = lnFactorialList[i-1].add(number);
		}
	}
	
	/*
	 * (non-Javadoc)
	 * @see statistics.AbstractStatisticTests#TestCasesOnly(int, int[][][], int[][][], int)
	 */
	public double TestCasesOnly ( int order, int[][][] observedCases ) 
	{		
		double z_cases = TestCasesOnly_Zscore(order, observedCases);
		double probability = Statistics.chiSquarePvalue(1, Math.pow(z_cases,2)) ; // / 2;
		return probability;
	}
	
	/*
	 * (non-Javadoc)
	 * @see statistics.AbstractStatisticTests#TestCasesVsControls(int, int[][][], int[][][], int, int[][][], int[][][], int)
	 */
	public double TestCasesVsControls ( int order, int[][][] observedCases, int[][][] observedControls ) 
	{ 
		double z_diff = TestCasesVsControls_Zscore(order, observedCases, observedControls);
		double probability = Statistics.chiSquarePvalue(1, Math.pow(z_diff,2)) ; // / 2;
		return probability;
	};
	
	
	public double TestCasesOnly_Zscore ( int order, int[][][] observedCases ) 
	{	
		int TOTAL_CASES = observedCases[1][1][1] + observedCases[1][1][0] + observedCases[1][0][1] + observedCases[1][0][0] ;
		double MIN_CASE_FREQ = (double) 1 / (TOTAL_CASES + 1);
		
		double P = (double) ( observedCases[1][1][1] + observedCases[1][1][0] ) / TOTAL_CASES;		
		double Q = (double) ( observedCases[1][1][1] + observedCases[1][0][1] ) / TOTAL_CASES;
		P = Math.max(P, MIN_CASE_FREQ);
		Q = Math.max(Q, MIN_CASE_FREQ);
		double expectedCoCarriers = P*Q*TOTAL_CASES;
		
		//No need to calculate the max : they all cancel out for case-only.
		double D_cases = (double) (observedCases[1][1][1] -  expectedCoCarriers) / TOTAL_CASES;		
		double sigma_cases  = Math.sqrt ( P*(1-P)*Q*(1-Q) / TOTAL_CASES );
		double z_cases = D_cases / sigma_cases;
		
		return z_cases;
	}
	
	
	public double TestCasesVsControls_Zscore ( int order, int[][][] observedCases, int[][][] observedControls ) 
	{ 
		int TOTAL_CASES = observedCases[1][1][1] + observedCases[1][1][0] + observedCases[1][0][1] + observedCases[1][0][0] ;
		int TOTAL_CONTROLS = observedControls[1][1][1] + observedControls[1][1][0] + observedControls[1][0][1] + observedControls[1][0][0] ;
		double MIN_CASE_FREQ = (double) 1 / (TOTAL_CASES + 1);
		double MIN_CONTROL_FREQ = (double) 1 / (TOTAL_CONTROLS + 1);
		
		//	Calculate case deviation
		double P = (double) ( observedCases[1][1][1] + observedCases[1][1][0] ) / TOTAL_CASES;		
		double Q = (double) ( observedCases[1][1][1] + observedCases[1][0][1] ) / TOTAL_CASES;
		P = Math.max(P, MIN_CASE_FREQ);
		Q = Math.max(Q, MIN_CASE_FREQ);
		double expectedCaseCarriers = P*Q*TOTAL_CASES;
		
		// r = D / sqrt ( P*Q*(1-P)*(1-Q) ) ~ Norm(0, 1/N)
		double D_cases = (double) (observedCases[1][1][1] -  expectedCaseCarriers) / TOTAL_CASES ;	
		double sigma_cases  = Math.sqrt( P*(1-P)*Q*(1-Q) ); 
		double r_cases = D_cases / sigma_cases ;

		//	Calculate control deviation
		double p = (double) ( observedControls[1][1][1] + observedControls[1][1][0] ) / TOTAL_CONTROLS;
		p = (double) ( observedControls[1][1][1] + observedControls[1][1][0] ) / TOTAL_CONTROLS;
		double q =  (double) ( observedControls[1][1][1] + observedControls[1][0][1] ) / TOTAL_CONTROLS;
		p = Math.max(p, MIN_CONTROL_FREQ);
		q = Math.max(q, MIN_CONTROL_FREQ);
		double expectedControlCarriers = p*q*TOTAL_CONTROLS;
		
		// r = D / sqrt ( p*q*(1-p)*(1-q) ) ~ Norm(0, 1/n)
		double D_controls = (double) (observedControls[1][1][1] -  expectedControlCarriers) / TOTAL_CONTROLS;
		double sigma_controls = Math.sqrt( p*(1-p)*q*(1-q) );		
		double r_controls = D_controls / sigma_controls ;
		
		//	correlation-contrast 
		double r_contrast = r_cases - r_controls ;
		double sigma_contrast = Math.sqrt( (double) 1 / TOTAL_CASES + (double) 1/TOTAL_CONTROLS ) ; 
		double z_contrast = r_contrast / sigma_contrast;	
		
		return z_contrast;
	}
	
}
